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Abstract 


The MOSTAS computer code for wind turbine analysis 
is reviewed, and the techniques and methods used in its 
analyses are described in some detail. Some impressions 
of its strengths and weakness, and some recommendations 
for its application, modification, and further develop- 
ment are made. Additionally, some basic techniques used 
in wind turbine stability and response analyses for systems 
with constant and periodic coefficients are reviewed in the 
Appendices . 
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Introduction 


The purpose of this memorandum is to briefly review the 
MOSTAS Code models and solution methods; to outline some im- 
pressions of its strengths and weaknesses; and to suggest 
some modifications and tests of the models and analysis methods 
employed. Also, this memorandum reviews some basic techniques 

used in wind turbine stability and response analyses. 

MOSTAS is a general computer analysis code for calculat- 
ing the dynamic loads, and for investigating the aeroelastic 

and mechanical stability of horizontal axis wind turbines and 
1 2 

helicopters. ' It was originally developed for helicopters 
and later adapted to wind turbines. This review will be con- 
cerned with its application to wind turbines. 

2:. Typical Horizontal Axis Wind Turbine Analysis 

To begin to evaluate MOSTAS, it is useful to imagine the 
steps one might take to perform a typical dynamic analysis of 
a horizontal axis wind turbine. 

1. Derive equations of motion (partial differential 
equations or finite element models) for blades, 
tower, pod, and other system components. Also, 
formulate aerodynamic, inertial, and gravity loads. 

2. Use equations of motion to calculate coupled normal 
modes (J)^, and natural frequencies, o)^. 


for the blades. 
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These would preferably be rotating modes, although 
non-rotating modes could also be used. Calculate 
also normal modes and frequencies for the tower 
and other components. 

3. Choose a sufficient number of modal coordinates 
to adequately describe the motions of the system. 

Then derive modal equations for the q^^ ' s from the 
equations of motion for the blades, tower, pod and 
other system components. 

4. Couple entire system together. 

5. Solve for the quiescent (static) blade response and 
corresponding static loads, bending moments, etc. 

6. Solve for the dynamic blade response to the time- 
varying aerodynamic, inertial, and gravity loads. 
Calculate the dynamic loads to be added to the static 
loads. Investigate for any dynamic instabilities, 
both aeroelastic and mechanical in origin. 

In attempting to carry out step 6, one often finds that the mass, 
stiffness, and damping terms for the blade equations con.tain 
periodic components in addition to the usual constant coefficients. 
For three or more bladed rotors, one can reduce most of these 

3 

periodic components by introducing multiblade coordinates. 

The equations can then be solved using standard, constant coef- 
ficient system response techniques. For a two or one bladed 
rotor, additional multiblade coordinates and harmonic balance 

methods need to be introduced to approximately eliminate the 
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. 4 

periodic coefficients. A crude method that has occasionally 

been used to deal with small periodic coefficients is to simply 
time-average the periodic variations over one blade revolution, 
thereby obtaining constant coefficient equations. However, un- 
less the periodic coefficients are small, this may result in 
some errors and may miss some potential areas of instability. 

As a practical alternative to introducing multiblade coordinate 
and harmonic balance methods, one may directly integrate the 
periodic coefficient equations numerically and introduce Floquet 
techniques to examine the system for instability.^'® 

The aerodynamic, inertial and gravity loadings on a wind 
turbine tend to occur periodically in multiples of the rotation 
frequency For constant coefficient equations, it is easiest 

to obtain the steady-state dynamic response to such loadings 
using frequency response techniques. These involve finding the 
harmonic response qj(t) = qje^^^m^ to each forcing frequency 

= mn , then summing up all the harmonic responses. See Appen- 
dix A. To investigate for instability, one recasts the dynamic 
equations in state-space form, determines the eigenvalues of the 
system, and notes if any eigenvalues are positive real, or have 
positive real parts. See Appendix A. For the direct numerical 
integration of the periodic coefficient equations, one picks a 
convenient integration scheme (Runge-Kutta, Newmark, Central 
Difference, etc.) and integrates from some initial condition. 

By proper choice of the initial condition, one can eliminate all 
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transients from the response and thereby show the desired steady- 
state dynamic response integrating through only one blade re- 
volution, instead of the very large number usually needed to 
reach steady-state for lightly damped systems. See Appendix B. 

To investigate these periodic coefficient equations for insta- 
bility, one uses Floquet techniques, which implies finding the 
eigenvalues of the Transition Matrix. See Appendix B. The use 
of multiblade coordinate and harmonic balance techniques is de- 
scribed in Appendix C. 

Having found the dynamic response by any of the methods 
outlined above, one may proceed to obtain the dynamic loads to 
be added to the static blade loads by summing all the additional 
applied aerodynamic, inertial, and gravity loads acting on the 
blade. This force summation method is preferable to the simpler 
mode deflection method since it represents the static loads 
more accurately, and it handles discontinuities better. Still, 
one should be careful to use a sufficient number of blade elastic 
modes to insure reasonable convergence in blade bending moments 
and shear. ^ 

The typical analysis of horizontal axis wind turbines out- 
lined here has several variations, depending on the authors. 

8 5 

See for example, Spera, Kottapalli and Friedmann, Warmbrodt 
6 9 

and Frxedmann, and Miller, et al. A good general review and 
bibliography for dynamic analysis of such horizontal axis wind 
turbines has recently been given by Friedmann. The specific 
steps and procedures used in MOSTAS will be described next. 
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3 . Description of MOSTAS Code 

3. 1 General Layout 

The MOSTAS Code is a very general computer formulation 
which attempts to include a wide variety of structural, inertial, 
aerodynamic, and generator system effects to obtain the dynamic 
loads on a wind turbine. The rotor may have one, two, or more 
blades which may be tapered, twisted, and preconed out of the 
plane of rotation. They may have cantilevered, hinged, or 
teetering attachments to the hub. The dynamics of the tower, 
pod, power train, and control system are considered in the analysis 

Basically, the MOSTAS Code can be divided into 3 sub-systems 
See Fig. 1. First, there is the MOSTAB-HFW system which calcu- 
lates the basic loads and bending moments on an isolated 
single blade assuming a fixed shaft (i.e., no pod or tower 
motion) , rotating at constant rotation speed Also, all the 

blades are combined to give total loads and moments f actinq 

O' 

at the fixed shaft. Additionally, time varying, linear math 
models for the isolated single blade are generated in MOSTAB-HFW. 
Second, there is the ROLIM system, which is used to assemble 
a linearized model of the rotor from the isolated single 
blade models generated in MOSTAB-HFW. This linear rotor model 
includes multiblade coordinate transformations and is used in 
subsequent coupled analyses of the wind turbine. Finally, 
there is the WINDLASS system which couples the ROLIM produced 
linear rotor model with linear models of the pod, the tower, 
the control system, and the power train. This WINDLASS system 
produces small perturbation loads f , which are to be added on to 
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the basic MOSTAB-HFW fixed shaft loads fj^ , to produce the 
final loads on the blades. The perturbation loads f contain 
all the effects of the shaft motions and the rotor interactions 
v;ith other components of the system. The input to this linear 
WINDLASS system are the fixed shaft, constant rotation speed rotor 
loads f^, previously found by MOSTAB-HFW. Thus, the overall 
logic of the MOSTAS code is to concentrate mainly on the fixed 
shaft loads, which probably comprise the greater part of the 
blade load, then to add on the generally small corrections due 
to shaft motions and subsequent pod, tower, control system, 
and power train system interactions. 

The above MOSTAS scheme here can be contrasted with some 

g 

of the other analyses quoted earlier. Spera, and Kottapalli 
and Friedmann^ essentially analyze the fixed shaft system and 
hence are comparable to the MOSTAB-HFW portion of MOSTAS. See 

g 

Fig. 1. Warmbrodt and Friedmann include the effects of tower 
and pod motions, and hence would be similar to the WINDLASS 
portion of MOSTAS, but with the steady wind, wind shear, tower 
shadow, gravity, etc. , input going directly into the rotor, 
rather than coming indirectly through the fixed shaft rotor loads 

9 

f onto the pod. Miller et al. consider the fixed shaft 
o 

loads by harmonic decomposition into static loads, and using 
approximate dynamic magnification factors to obtain the dynamic 
loads. The effects of the tower motions are found by analyzing 
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simple tower-pod- rotor models. 

The MOSTAS treatment of teetering rotors should also be 
mentioned. In order to calculate the fixed shaft loads fj^ 
for such a teetering rotor, the teetering angle y of the' rotor 
must be known. Since this y depends on the interactive motion 
of all the blades, a separate linearized analysis of the rotor 
blades is required. This is done as an additional small sub- 
loop in the MOSTAB-HFW system shown in Fig. 1. A linearized 
model of the rotor blades with teeter is constructed within 
MOSTAB-HFW; these are solved by a coupled linear analysis for 
the linear blade deflections and teetering angle; then this 
teetering angle is used to modify the total rotor load f^ 
feeding into the WINDLASS portion of the analysis. Fig. 1. 

This procedure seems a complex way to include the additional 
teetering degree of freedom y, but it does adapt to the frame- 
work of the isolated blade model scheme of MOSTAB-HFW. 

More specific details of the various MOSTAS component sys- 
tems are described next. 

3.2 MOSTAB-HFW System 

The MOSTAB-HFW system deals with obtaining the dynamic 
loads and response of an isolated blade assuming a fixed shaft 
(i.e., no pod or tower motion) , rotating at constant rotation 
speed It is essentially a modal analysis, and the blade 

motion is described by a superposition of blade natural vibration 
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modes added to a quiescent (static) blade position. 

MOSTAB-HFW uses as input one to four blade natural modes 
and frequencies, the blade mass characteristics, the aerodynamic 
section characteristics, and the quiescent shape of the blade 
(the static shape before dynamic vibrations occur) . It first 
proceeds to characterize the quiescent shape of a blade reference 
line (BRL) and a set of blade axes in terms of 3 translations and 
3 Euler rotations from a set of rotor reference coordinates which 
rotate with the rotor. This is carried along numerically by a 
6x1 column vector w^(s) for every point s on the BRL. Then 
it considers the blade vibration modes as small perturbations 


Aj(s) 3 j (t) about the quiescent position w^(s). Each mode shape 
Aj(s) is again carried along as a 6 x 1 column vector, and by 
assuming the orthogonality properties of the given modes, the 


modal equations of motion can be expressed in the standard form, 

ft' Pi = j A- ^ is (1) 

o 

where, -p 

c/s 

— Blade natural frequency 


In the above procedure, all geometric nonlinearities are retained 
in characterizing the quiescent position of the blade w^(s), but 
small perturbation deflections are assumed for the subsequent 
vibrations. The given mode shape data Aj(s) and natural fre- 
quencies (jjj are usually taken as rotating blade modes (centrifu- 
gal force effects included) . The external perturbation force p 
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consists of all external perturbation loads except for the 
inertia loads and centrifugal stiffening loads (these have 
already been accounted for on the left-hand-side of Eq. (1) ) . 

Rather than evaluating the generalized mass and generali- 
zed force integrals in Eq. (1) directly, the MOSTAB-HFW pro- 
gram prefers to mechanize all computations around repeated 
evaluations of the following integral, 

$ = I a’[cs) P(pJ , p, s, t) ds <2) 

o 

where S is a desired integral and P is the total load on the 
blade from all sources (aerodynamic, inertial, centrifugal, 
gravity, etc.), and includes both linear and nonlinear effects. 
The external perturbation load p. in Eq. (1) is related to the 
total load P by the relation, 

^ =1 P ~ (2-^ + p3 + Aj (3) 

where p^ is the quiescent (static) load on the blade, while 
mAj3j and are the inertia and centrifugal stiffening 

loads that must be added in order to cancel out the corres- 
ponding loads -hA-. B- and A\ already present in the 

total load P. 

The generalized mass term in Eq. (1) is obtained from 
the general integral formula Eq. (2) by numerical differentiation. 
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First Eq. (2) is evaluated with all motions 3 = 0 = B = 0 
and all external perturbation loads p = 0 in order 
to obtain the value, ■=■ J A- l^o • (See Eq. (3) 

rearranged to solve for P) . A subsequent evaluation with 

m m 

(3 • •= € (where £ is a small positive value) and 

all other 6 = ^ = P =0 will yield ~ 6 j Aj Aj ds 

Then , 

M;; = - (4) 

e 

A similar numerical dif ferentiation^ procedure with Pj ~ ^ 

» 0 * 

and all other ^ ^ = ^=0 will yield the linear stiffen- 
ing terms . Finally, the generalized force in the 

basic Eq. (1) is evaluated by placing the p given by Eq, (3) 
into the integral formula Eq.(2) to obtain the right-hand-side 
of Eq. (1) as, 

. M'‘ JaJp- ds = m''[S-^.] + § + m''Ki§ 

o 

In the above, all matrix exoressions M, S, S , K,. have been 

o I 

obtained by using repeated evaluations of the general integral 
formula Eq. (2) , The 3 term appearing in Eq. (5) is present in 
order to cancel the corresponding inertia forces present in 
the total load term S. 

The aerodynamic loads in the total load P are calculated 
assuming quasi-steady strip theory, and static stall effects are 
included by using curve-fit functions, ^^^^us angle 
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of attack. Geometric nonlinearities are retained since angles 
of attack are found numerically. The inflow of air may vary 
over the disk in accordance with simple momentum theory which 
treats independent annular volumes of air. Also, inflow cor- 
rections can be made for wind shear, tower shadow, and rotor tip 
losses . 

The equations of motion, Eq. (1) with the right-hand-side 
expressed by Eq. (5) are of the general form, 

p. + W p.; = ( Pj , pj, §: ,t) (6) 

These equations are integrated numerically around one rotor 
revolution, from 'f' — =0 to ^ — 2.TT , for a given 

set of initial conditions on ^ and ^ . The correct initial 

conditions that will result in steady-state , periodic solutions , 
are found using the procedure described in Appendix B. The 
numerical integration method used in MOSTAB-HFW consists of 
assuming g^^ to be constant over some small azimuth interval 
A'f' = 5tA‘t . The equations are then uncoupled, and a conven- 
tional solution over that interval is given as. 



( 7 ) 
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Care must be taken that the time interval Atr — is suf- 

ficiently short for the highest natural frequency CO,, involved, 
otherwise numerical instability results ( Air < ) • 

Some schemes for suppressing the numerical instability are dis- 
cussed, so that excessively small azimuth intervals (and hence 
excessive computer time) need not always be used. 

• • • 

Having obtained the required time histories ^ over 

one rotor revolution, the desired shears and bending moments f, 

b 

for an isolated blade are found from the relation, 

a 

P( P, p/p, ^ ,t) ,7) 

where f^ is a 6 x 1 column vector representing the 3 total 
shears and 3 total bending moments at any point s on the blade, 
and P is the total load on the blade from all sources (aerody- 
namic, inertial, centrifugal, gravity, etc.). The above rela- 
tion Eq.(7) is essentially a "force summation" method instead 
of a "mode displacement" method for the loads, and it is evalu- 
ated using the familiar set-up for the general integral formula 
Eq. (2) , only now Aj^(s) is replaced by some other function 
R (s^Tl") f and the lower limit o is replaced by the more general 
value s. 

The isolated blade load time histories f^^ computed above 
are additionally fourier-analyzed to obtain the amplitudes and 
phase angles of the lowest few harmonics involved. The loads 
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frcn'. each individual blade can then be combined together to 
form the fixed shaft load input f^ to the WINDLASS system as 
described in Section 3.1. 

Additionally, time varying, linear math models for the 
isolated single blade are generated in MOSTAB-HFW. This is 
accomplished apparently by combining the modal equations of 
motion Eq. (1) , with numerical differentiation of the general 
integral formula Eq. (2) , to obtain *3/ 3/ and 3 coefficients 
similar to the manner described in Section 3.2 for obtaining 
the generalized mass term This is done at many azimuthal 

locations of the blade. Because of the numerical differentiation 
all nonlinear effects are included in the derivative calculations 

Finally, it should be mentioned that for the case of a 
teetering rotor, the calculation of the fixed shaft loads f^ 
is considerably more complicated since the interactive motion 
of all the blades must be considered. This involves doing, 

(in the MOSTAB-HFW system) , a coupled linear analysis of the 
teetering degree of freedom y with the given modes of all the 
blades, as was described briefly at the end of the previous 
Section 3.1. 

3.3 ROLIM System 

From details inferred in Refs. 1 and 2, the ROLIM system 
basically assembles a linear (perturbation) math model of the 
entire rotor about the trim operating condition, from the linear 
math models of the isolated single blades generated in the 
MGSTAB-ilFW analysis. 
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The isolated single blade model is expanded to 
represent the full N-bladed rotor by assembling each set 
of modal coordinates in a stacked fashion. 


r 



\ 


if- 

I i 


where are the M modal coordinates of the blade, 

and the coordinates of each blade are assumed to be shifted 

in azimuthal phase an amount Aij; = 27r/N from its neighbor. 

(k) 

In addition to the above blade coordinates 3^ , one has 6 

additional rotor shaft coordinates , representing possible 
displacements of the rotor hub in three directions, and pos- 
sible rotations of the rotor hub about three axes due to 
motion of the tower and pod . The equations of motion are ex- 
panded to include these additional hub translation and rotation 
effects. Because the rotor hub motions are described in a 
fixed reference frame while the blade motions are described 
relative to a rotating reference frame, the resulting equations 
of motion will generally have mass, damping, and stiffness co- 
efficients which are now periodic functions of azimuth position 
of the k^^ blade, 
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For rotors with 3 or more blades, it is known that the 
periodic coefficients in the above equations can be largely 
cut down and reduced to constant coefficients by introdu- 


cing multiblade coordinates 


isi' "lei' 2si' 


such that, 


p*- p, 






where = fit + (k-1) 2 it/N represents the azimuthal position 
of the blade and the number of coordinates taken equals 

the number of blades, N. See Appendix C. The ROLIM system 
introduces such multiblade coordinates and resulting equations 

of the entire rotor are recast in the matrix form. 


Mu. + Q ^ c, j S2j 


where M, P, Q, are constant coefficient square matrices, 
and y is a column matrix of the unknown coordinates. 


i ' J 


The first 6 coordinates x^ in y above represent the rotor shaft 
motions while the remaining b^'s, etc. are the multiblade 
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coordinates. The right-hand-side R represents the effects 
of forces and moments acting on the rotor and is considered 
a function of, 

f - forces and moments from nacelle pod to rotor hub. 
c - control inputs to the blade angles. 

- speed of rotation, 
g - wind velocity vector. 

The above linear model of the rotor system Eq. (9) , is then 
used in the subsequent rotor-tower interaction system 
WINDLASS described in the next section. 

The second order Eqs. (9) can also be rewritten in the 
state-space form as twice as many first order equations, 

P - Q ^ ^ 


where the matrices P, Q, R, and y now have different defi- 
nitions and dimensions. See Appendix A. ROLIM assembles some 
of its analyses in the form of Eqs. (11) as well as Eqs. (9). 

Some comments should be made about the use of multi- 
blade coordinates when the number of blades N < 3. For 
2-bladed rotors, N =■ 2 , the multiblade coordinates Eq. (8) 
will not eliminate the periodic terms. Rather, it will 

introduce the two coordinates b_ . (t) and b, . (t) which toaether 

Ti Ai 

with x^ , can then be expanded in a harmonic series, 
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where x . 

01 

' ^si' 

^Toi' ^TlSi' •• 

. are all functions of 

time . ■ 


The harmonic balance method can then be applied to obtain 
approximate solutions as discussed in Appendix C. The re- 
sulting response will then consist of a fundamental frequency 
plus various harmonics of frequency m^/ where m is an integer. 

The ROLIM system introduces a special form, for the multi- 
blade coordinates for a two-bladed rotor, namely, 

^ • . / li 


where for N = 2, one has + (k-Dir. In matrix form 

this appears as. 




0 



i 

Slcc^^^ 


0 



1 





(14) 


This is placed into Eqs. (6) written in state-space form (i.e., 
as first order equations), and after adding the hub motion 
coordinates and multiplying by the transpose of the above 
transformation, one obtains equations of the form of Eq. (11) 
where y now represents the variables x . , k. , b , b , b , b . 
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The transformation Eq. (14) is a mathematically valid one, and 

it appears efficient since it represents the 4 unknowns 6^^^, 

^(1)^ g(2)^ ^(2) terms of 4 new unknowns b^, b^, b^, b^, and 

also, it tends to change the sin variations into sin 

2 

variations. However, Kaza, Janetzke and Sullivan pointed out 
by •opolioation to a simple example, that if one subsequently 
time-averages out the . resulting sin and cos 2\p^ variations 

to zero, it gives poor results for the frequencies and instability 
regions . 

A somewhat better transformation to use for these 2-bladed 
rotors is the following. 


or in matrix form, this appears as. 




1 

0 


0 


0 



0 

I 





1 


1 

0 


0 


o 

L P J 

i 

0 

1 



SI 







( 16 ) 






I 

This transformation represents the 4 unknowns 6^^^, 

• ( 2 ) ... 

3 in terms of 6 new unknowns b,b,b,b,b,b, and 

o' o s s c c 

upon substitution into the equations of motion and multiplying 
by its transpose, it will lead to 6 blade equations instead of 
the 4 obtained by the previous transformation, Eq. (14). Although 
this does not appear as efficient as the previous transformation. 
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it allows more flexibility among the coordinates and gives much 
better results if one subsequently time-averages the resulting 
sin and cos 2Tj;^ variations to zero. It gives reasonable 
frequencies and it will preserve the instability regions 
(slightly shifted). See Refs. 4 and 14. 

It would appear then, that ROLIM would do better using 
the transformation Eq. (16) with time-averaging to obtain constant 
coefficient equations rather than using the special form, Eq. (14) 
with time-averaging. Otherwise, if the special transformation 
Eq. (14) is used, one should solve the resulting equations and 
check their stability using Floquet techniques, as discussed 
in Appendix B, or by some other scheme which includes the 
periodic coefficients in the equations. This was also suggested 
in Ref. 2 . 
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3»4 WINDLASS System 


The WINDLASS system couples together linear models of 
the rotor, the nacelle pod, the tower, the control system, 
and the power train, to give the dynamic response of the 
overall system. The main input to the WINDLASS system is 
the fixed shaft, constant rotation speed loads f^, previous- 
ly found by MOSTAB-HFW and applied to the nacelle pod. This 
sets up small perturbation responses everywhere in the over- 
all coupled dynamic system. The total loads on the blades 
are then taken as the sum of the fixed shaft loads f, and 

D 

the additional perturbation loads f due to the nacelle pod 
motions. Figure 2, taken from Ref. 1, shows a block diagram 
of the overall WINDLASS system, its five subcomponents, its 
variables, and the general form of its equations. A brief 
description of the various subcomponents follows. 

The linear rotor model is developed by the ROLIM system 
from the MOSTAB-HFW analysis and was described in the previous 
section. The governing equations are given by Eqs. (9) or (11) 
and the variables y are given by Eq.(lO). The first six of 
the variables y represent the six absolute rotor shaft motions 
X, while the remaining variables represent the multi blade co- 


ordinates b.-. , b . , . . 

I iL S IL 

inputs to the equations 


for the rotor blades. The various 
f, c, (7, g are also described in the 


previous section on ROLIM. 
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The nacelle pod acts as an interfacing device that con- 
nects the rotor, tower, control system, and power train sys- 
tems together. The pod itself is modeled as a massless elas- 
tic spring superimposed with an infinitely rigid mass, so that 
the pod model has no relative structural vibration modes, but 
it does contribute its mass properties to the overall system 
dynamics. The yaw drive spring stiffness plays an important 
role in the nacelle pod representation. The form of the govern- 
ing equations for the pod can be seen in Fig. 2. The fixed shaft 
rotor loads f^^ are the primary forcing function on the pod, while 
f is the perturbation load from the pod to the rotor. The f^ 
are loads from the pod to the tower top, x^ is the absolute 
position of the tower top, f^^ are the loads applied to the pod 
rigid body mass, x^^ is the absolute position of the pod rigid 
mass, and Yd is the torque from the rotor to the power train. 

The remaining subscripted f, x, and y variables represent various 
intermediate loads, positions, and torques. 

The tower model is represented by its natural modes of 
vibration These are found from a separate finite element 

analysis of the tower, and are usually performed using a mass 
at the top to approximate the mass properties of the nacelle- 
rotor unit. The resulting modes and frequencies then’ give a more 
accurate description of what the tower is doing during its 
system coupled motion. However, the effect of this mass has 
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to be subtracted out from the model of the isolated tower used 
in the WINDLASS system. Also, the effect of the elastic tower 
mounted on a flexible base is also included so that one might 
additionally consider a movable soil. Six additional degrees 
of freedom 6, representing rigid body motions of the base, are 
included in the tower model, and appropriate equations are de- 
veloped. The form of the governing equations for the tower 
are shown in Fig. 2, one set representing the predominantly 
elastic tower vibration modes 5, the other representing the 
predominantly base vibration modes 6. The f^ are the loads 
from the nacelle pod to the tower top, the is the absolute 
position of the tower top, and the f„_ are external applied 
loads to the tower such as due to drag or cables. It would 
seem the tower model could probably be simplified by includ- 
ing flexible soil effects as a simple reduction in the tower 
natural frequency, instead of introducing six additional degrees 
of freedom. 

The control system model represents the power machinery, 
power machinery controls, utility network dynamics, rotor speed 
controller, and other servo-systems present to regulate the 
speed of the wind turbine. A simple such model is shown in 
Fig. 2. Here h represents the command inputs to the system; 
i.e., the desired speed or the speed governor, the a represents 
the control system degrees of freedom, and 4> and <(> represent 
the angular position and velocity of some gear of the rotor shaft. 
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The controller generates a rotor control input c from the 
control system degrees of freedom a, which is then fed into 
the rotor to control its behavior. The control input c 
usually represents a rotor collective pitch change 0^. The 
associated torques on the power train and nacelle pod 
are also generated from the control system degrees of 
freedom a. 

The power train model represents the assembly of gears, 
shafts, and generator that make up rotating parts of the wind 
turbine. The WINDLASS system makes up the power train by 
connecting a series of basic modules, each of which consists 
of a large main gear, a small pinion gear and an interconnect- 
ing elastic shaft. The smaller pinion gear is considered mass- 
less, and two damping coefficients are considered, one for the 
shaft and one for the bearings. Additionally, the torques 
transmitted through the system and to the supporting bearings 
are also considered. The form of the governing equations for 
the power train are shown in Fig. 2. Here (ji represents the 
power train degrees of freedom which relate to the angular 
position of each gear, while y„ represents the sum of the torque 
acting on the power train coming from the rotor from the con 

trol system y^, and from externally applied torques Ygp- The 
rotor speed is related to ^ , and the y represents the base 

ID 

torque loads on the pod. 
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The equations of motion for the five subcomponents of the 
WINDLASS system and their interrelationship is shown in Fig. 2. 

It remains to assemble all the equations together into a con- 
venient form to solve them simultaneously. Since the wind 
turbine system was developed, for the sake of generality, as a 
series of independent blocks linked together, the equations as 
presented include displacements and interface force terms as 
variables. To solve such a large set of equations is unwieldy 
since not only are there many variables, but also the deflec- 
tion and force operators may contain terms of widely varying 
magnitude which may cause numerical problems. To reduce the 
size of these equations, the WINDLASS procedure separates the 
variables into two groups - an "independent" group, usually 
comprising the displacements, and an "eliminative" group, usually 
comprising the interface forces. The equations are then also 
separated into independent and eliminative equations so that 
they can be written as. 


M ^ W + W + Q. + V 

Eg e = Eyv ^ + E;; w + V 


( 17 ) 


where w represents the independent variables, e represents the 
eliminative variables, and v represents the externally applied 
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loads. The e variables can then be eliminated by solving the 
second equation for e and placing into the first equation to 
yield a smaller set of equations, 



3 W 
9 " 



(18) 


where , 





S 


, etc. 


After solving the smaller set of Eqs. (18) for w, the elimina- 
tive variables e are obtained from the second equation Eq. ( 17 ) by 
multiplying by E^^ . It should be noted that in the above 
reduction procedure, the eliminative variables e never appear 
as time derivatives in the equations of motion. Hence the 
order of the differential equation system in Eqs» (17) is the 
same as that in Eq. (18) even though Eq- (18) has fewer vari- 
ables. Thus no essential dynamics has been lost. Also, it 
should be noted that this reduction procedure requires the 
matrix E^ to be nonsingular. Often, because of the choice of 

variables and the existence of constraints, E^ is singular. In 
such cases, a special processing technique is used to get around 

this by first using linear, algebra techniques to express the 
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raatrix E in the form, 
e 

Ee = C"' D B 


( 19 ) 


where C and B are nonsingular matrices and D is the Identity 
matrix except for zeros in the first n rov;s corresponding to 
the order of the singularity in E^. This form isolates the 
singularity and allows one to rewrite Eqs. (17) in the same 
form as before, only now w has n additional variables, the 


M,B,K,V,E,E*,E.. ,E ^ 

g' g g' g' w' w' w v matrices are defined somewhat 

differently, and E^ becomes the unit matrix. The reduction 

to the smaller set of Eqs. (18) then follows. 

The WINDLASS system can solve the assembled linear equa- 
tions of motion Eqs. (18) by either using frequency response 
methods or by using direct numerical time integration methods. 
Originally, only the frequency response method was available, 
but later, the time integration method was added. 

In the frequency response method, the forcing functions v 
v/hich generally occur periodically over a cycle of turbine 
revolution, are broken down into their harmonic components. 

The harmonic responses w are then found and added together to 
give the total response, as shown in Appendix A. The WINDLASS 
subprogram to obtain the frequency response, Eqs. (A- 8) and (A-9) 
is designated DYNAM2 , while the subprogram to recover the time 
history Eq. (A-10) is designated REC0V2. The total loads acting 
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on the blades fipQip a^re then found by summing the original 

fixed loads fj^ and the here . calculated perturbation loads 

f. Also if desired, the blade response relative to the ro- 

(k) 

tating system B can be obtained from the multiblade co- 
ordinate transformation Eq. (8). It should be noted that 
this frequency response method as discussed in Appendix A 
requires that the Eqs. (18) have constant coefficients. If 
the equations have periodic coefficients, the periodic parts 
must be time-averaged over one cycle, and only the constant 
portions used. This may lead to some errors for 2-bladed 
rotors, particularly if the specialized 2-bladed multiblade 
coordinate transformation, Eq. (14) is used. See Ref. 2. 

The overall MOSTAS system using DYNAM2 and REC0V2 is labeled 
as MOSTAS-A in Ref. 2. 

The time integration method of solving Eqs. (18; is con- 
venient for arbitrary time varying loads applied to the struc- 
ture (such as wind gusts) , and it allows any periodic coef- 
ficients in Eqs. (18) to be accounted for. The equations of 
motion Eqs. (18) could be numerically integrated by any of 
the standard numerical techniques such as Runge-Kutta, Newmark, 
Central Differences, Houbolt, etc., to get the response. How- 
ever, because of the very general nature of the assembly methods, 
there may be constraints among the various chosen coordinates w 
such that they are not truly independent, and hence the mass 
matrix M may be singular. To get around this singularity, 
WINDLASS makes use of eigenvalue analysis to arrive at the 
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following numerical integration method. First, it rewrites 
Eqs. (18) in first order form, 

T" Cy — Q ly Rv (20) 

where the matrices P, Q, R, v, y, are defined by comparing 
with Eqs. (A- 2) and (A-1) of Appendix A. Then, it separates 
out the time-averaged constant parts P^, from the periodic 
parts P, Q, of the P, Q matrices and rewrites Eqs. (20) as 

- Kv - F Cy + ( 21 ) 

where the entire right hand side will be treated as a constant 
forcing function over small time intervals At in the subsequent 
numerical integration process. ' Next, using eigenvalue analysis, 
the program finds the eigenvalues and eigenvectors y^, 
of the homogeneous system P^y -Q^y = 0 and of its transpose 

T • ■ T 

P z -Q z = 0.: Sxnce the matrices P and Q may in general 

o o i o o ^ 

be singular, the eigenvalues and eigenvectors y^ and are 
found using special linear algebra methods described in Ref. 11 
which first take out the singularities in P^ and Q^. The re- 
sulting eigenvectors y^ are arranged, column by column, in an 
array Y called the modal matrix. Similarly, the z^ are arranged 

into another modal matrix Z . Then introducing new coordinates 

T 

q such that y = Yq and multiplying Eq. (21) by Z , the program 
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where and is generally complex. This numerical 

integration scheme is the first order complex counterpart of 
that used in MOSTAB-HFW involving second order real equations, 
Eqs. (6) and (7) . 

Finally, after solving for q, the physical coordinates 
are obtained through the transformation y = Yq, and the time 
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histories of the blade loads f.pQ.p are obtained as before by 
summing the original fixed shaft loads fj^ and the here calculated 
perturbation loads f. In the above procedure for integrating 
Eqs. (18) / it should be mentioned that for practical expediency, 
the number of degrees of freedom is substantially reduced prior 
to the solution. This is done by discarding the very high fre- 
quency modes in Eqs. (18) . Otherwise, time integrations of these 
equations would be extremely difficult. 

The WINDLASS subprogram that performs the numerical time 
integration method here is designated WINDGUST, and can be used 
to obtain the response to any forcing function v in Eqs. (20) 
subject to any initial conditions. For forcing functions which 
occur periodically over one cycle of turbine revolution, the 
resulting steady-state response can be found by integrating over 
only one cycle rather than many cycles provided the proper initial 
conditions are found. WINDGUST uses the scheme described in 
Appendix B for this purpose. The overall MOSTAS system using 
the time integration WINDGUST subprogram is labeled as MOSTAS-B. 

Summarizing the two methods of solving the linear equations 

of motion, MOSTAS-A uses frequency response methods (DYNAM2 and 

RECOV2) to solve Eqs. (18) with time-averaged coefficients, while 

MOSTAS-B uses time integration methods (WINDGUST) to solve Eqs. 

(21) with the periodic coefficients P and Q. From an investi- 

2 

gation by Kaza, Janetzke, and Sullivan , MOSTAS-B seems to do 
better for the forced response than MOSTAS-A when applied to 
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2-bladed rotors, since it takes into account the periodic co- 
efficients P and Q of Eqs. (21) . On the other hand, MOSTAS-A 
seems a simple and efficient method of getting the harmonic 
response when periodic coefficients are small, as in 3-bladed 
turbines. 

For investigating the stability of the rotor systems 
given by the assembled linear Eqs. (18) , the WINDLASS system 
sets the forcing functions v = 0, then time-averages the co- 
efficients to constant values, then obtains the eigenvalues 
p^ of the system and notes if any are positive or have posi- 
tive real parts. See Appendix A. Because the mass matrix M 
is often singular, the methods previously described for the 
WINDGUST analysis, Eqs. (20^ to (23) are used to get around 

this singularity. In fact, the eigenvalues 

found from Eqs. (22) are exactly the desired eigenvalues p^^ 
to be examined for stability. Hence, the previous WINDGUST 
response analysis automatically includes the stability in- 
vestigation of the time-averaged constant coefficient system. 

To examine the stability of Eqs. (18) when periodic co- 
efficients are present, Floquet methods should be introduced 
as described in Appendix B. It does not appear in 
Ref. 1 that this has been incorporated into the WINDLASS 
analysis , but in any event, its inclusion is 

simple since it involves merely finding the eigenvalues Aj, 
of the "transition matrix" Q defined in Eq. (B-5). This 
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matrix Q has already been formed by WINDGUST to obtain the 
proper initial conditions for steady-state response. Appen- 
dix B. If any of its eigenvalues has a magnitude equal 
to or greater than unity, the system is unstable. Alterna- 
tively, a root perturbation method is currently being considered 
for investigating the stability of the periodic coefficient 
systems. 

As mentioned earlier, the use of the time-averaged 
constant coefficient method should be used with caution when 
dealing with 2-bladed rotors. See Kaza, Janetzke, and Sullivan", 
and the discussion earlier in Section 3.3. This is apparently 
due to the nature of the 2-bladed transformation, Eq. (14) 


used in ROLIM. 
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4 . Strengths and Weaknesses of MOSTAS 

The MOSTAS code for wind turbine analysis seems to be 
a very general code incorporating some interesting features, 
which has been built up over the years to handle a variety 
of helicopter and wind turbine problems. Many additions 
have been made to the original base to make the routine more 
general, and many sophisticated data processing technqiues 
have been added accordingly to handle these. As a result, 
one begins to lose the feel of what is going on in the struc- 
ture, and one wonders if numerical accuracy is beginning to 
obscure even simple results. 

The following are some specific comments about its 
strengths and weaknesses as they appear to the authors of 
this review. 

Strengths 

1. Obtains blade loads by finding basic isolated blade 
loads plus smaller linear perturbation loads due to 
tower motions. 

2. Works with complete loads on blades from all sources, 
(aerodynamic, inertial, centrifugal, gravity, etc.) 
Can include stall, variable inflow, and other aero- 
dynamic nonlinearities. 


3 . 
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4. Modular arrangement of components. Can handle all 
kinds of subpieces. 

5. Good coupling and assembly of systems together. 

Uses good data processing methods. Numerical in- 
tegration methods seem reasonable. 

6 . Includes frequency response methods for harmonic 
loads as well as time integration methods for load 
time histories. 

7. Good procedure for obtaining correct initial condi- 
tions for steady-state solution. 

8. Includes effects of periodic coefficients for forced 
response and possibly for stability investigations. 

VJeaknesses 

1. Seems a reasonable, but somewhat cumbersome program. 

Very long, does everything numerically, most struc- 

* 

tural nonlinearities retained, numerical differen- 
tiation may lead to small differences of large num- 
bers, difficult to see what is causing what. 

2. Many small couplings and effects are included which 
are probably unnecessary, but numerically complica- 
ting. Because allowance is made for so many features, 


•k 

Wind turbine blades are generally stiffer than helicopter 
blades. Can neglect many of the nonlinearities retained 
here by order of magnitude assumptions. 
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the program does things very generally and requires 
much data processing and sophisticated techniques 
to get answers . 

3. Method is essentially a modal analysis for the blades. 
Use of at least 4 blade modes may be required for 
finding bending moments and shears on blades. This 
may tax the numerical capability of this program. 

4. Quiescent position not consistently found in the 
analysis. Must get from elsewhere. Seems one is 
using very sophisticated additions to a crude model. 

5. Need separate programs to obtain modes and quies- 
cent position for the analysis. It seems the code 
should be able to get these for you, somewhere. 

6. Teetering system seems inefficient in the operation 
of the system. Could perhaps be improved. 

7. Some problems exist with the use of multiblade coordi- 
nates for 2-bladed rotors. The special form used in 
the program may be changed as discussed. 

8. Possibility of flutter instability when there is a low 
torsional flexibility of the blade due to a loose 
pitch-change mechanism. Should be investigated more. 
Some question whether the system can model blade aero- 
elastic instabilities properly, since unsteady aero- 
dynamic forces and wake interference effects may have 


to be included. 
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5. Recommendations 

This report has attempted to give a review of the MOSTAS 
computer code for wind turbines, and to describe some of the 
techniques and methods used in its analyses. Also, some 
general methods used in wind turbine stability and response 
analyses have been reviewed. Based on this study, the follow- 
ing recommendations are suggested. 

1. Make smaller, simpler models involving specific sub- 
components, to investigate the main origin of large 
loads and instabilities. Try to understand what 
causes what by these simpler models. 

2. Check out MOSTAS versus some simpler models to see 
what can be eliminated and simplified. 

3. Look more into effects of aeroelastic flutter insta- 
bilities and also mechanical instabilities on these 
systems, especially with proposed soft, flexible moun- 
tings . 

4. Look more closely at teetering effects and propeller 
whirl type flutter of these systems. This involves 
aeroelastic coupling with the tower yaw and pitch 
mechanisms . 

Change the special form of the 2-bladed multiblade 
coordinate transformation used in the analysis, or 
use Floquet methods (or their equivalent) to analyze 
the stability of the response. 


5 . 
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6. Improve teetering set-up if possible. 

7. Look at drive train and torsional shaft-generator 
coupling problems more closely. These can probably 
be uncoupled from much of the tower motion dynamics. 
Should calculate the quiescent position by the 
computer code, for the given geometrical blade shape. 


8 . 



Appendix A 

Stability and Response of Constant Coefficient Systems 


Given a system of N linear differential equations with 
constant coefficients, 


where M, B, and K are the square matrices of order NxN, while 
and F(t) are column matrices of order Nxl. These can be re- 
arranged as, 



Then, multiplying through by the inverse of the mass matrix 
gives 2N first order equations. 







(A-3) 


where A is a square matrix of order 2Nx2N, while and G are 
column matrices of order 2Nxl given by 



0 

-M*' K 





The above rearrangement, Eq. (A-3) , is valid providing the mass 
M is not singular, which is usually the case with physical 
systems. 

(a) Stability 


To investigate stability, one sets F = 0 (which gives 
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G = 0) in Eq. (A-3) to obtain a set of homogeneous equations. 
Then one seeks exponential solutions of the form v. ~ Y. 
Placing these into Eq. (A-3) leads to the standard eigenvalue 


problem, 


A ^ 


(A-5) 


Eigenvalues Pj^ of the matrix A can be obtained by standard 
numerical eigenvalue routines. If any eigenvalue Pj^ is positive 
real or has a positive real part, the system represented by 
Eq. (A-3) or equivalently by Eq. (A-1) is unstable. 

(b) Forced Response 

Under steady-state conditions, the forces F(t) on a rotating 
system tend to occur periodically in multiples of the rotation 
frequency (2. One can then express the force for a particular 
frequency = mfi, in the form, 

F (t) — {Ra. ^ F e O^t — 

(A-6) 

The response 3 ^(t) is similarly of the form, 

^(O = (A-7) 

Placing Eqs.(A-6) and (A-7) into the basic Eq. (A-1) and matching 
sine and cosine terms gives a set of 2Nx2N real equations. 


ti a 1 


(A-8) 


where one has the matrix elements. 



-40- 


£ ::r K ~ M 


, H = B 


(A-9) 


Given the amount of the m^^ harmonic force present F and 

, Eq. (A-8) can be solved by simple inversion to find the 
response arid each harmonic. Then, one may sum 
up all the harmonics to give the total periodic response as, 


C^cr^ 




Finding the response c[(t) this way rather than by direct numerical 
integration, allows one to assess the effects of a particular 
harmonic on the resulting response of the system. 
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Floquet Methods for Periodic Coefficient Systems 


Assume the coefficients M, B, K in Eq. (A-1) or equivalently 
the coefficients A in Eq. (A-3) vary periodically in time, rather 
than being constants. For illustrating Floquet methods, it 
will be convenient to use the first order representation, namely 
2N equations of the form, 

^ ' Aw ~ (B-1) 


where A(t) and G(t) are periodic over an interval T. 

(a) Stability 

The Floquet stability analysis described here follows that 

12 

given by Peters and Hohenemser . To investigate stability, one 
sets G=0 in Eq. (B-1) to obtain homogeneous equations. The Floquet 


theorem states the solution of Eq. (B-1) with G=0 is of the form 


= 3(t) 


(B-2) 


r k 1 

’-’here y_(t) and } are 2Nxl column matrices, and B(t) is 

' 2Nx2N square matrix periodic over period T, that is, B(T)=B(0) 


From the above, one can express 


x(°) = B (o) ^ 

C 7 C 


(B-4) 
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Also, one can express ;^{T) asj 




(B-s: 


II 

[q1 


(1) , 


where y . is the solution at t=T of Eq. (B-1) with G=0, for 

12 ) 

the initial conditions ^il remaining y^(0)=0, y 

is the solution for y 2 ( 0 )=l and all remaining y^(0)=0, etc. 
The square matrix [Q] is called the "Transition Matrix. " 
Equating Eq. (B-5) to (B-4) and introducing Eq. (B-3) gives, 

[c?][[bw} c, * 


C, -t - 


= [b(o)| c, &' i- 




(B-6) 


Fr' 

e. + 


Since Cj^ are independent, one must have 




P, T 
k 


(B-7) 


where the eigenvalues of the [Q] matrix. One then 

has the relation 




(B-8) 


from which the real and imaginary parts of the stability exponent 
are given as 


(B-9) 
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— Ip (B-10) 

The real part is a measure of the growth or decay of the 
response, as can be seen from Eq. (B-2) . Values of otj^>0 (or 
equivalently |Xj^|>l) indicate instability. The imaginary part 
represents the frequency. However, because tan~^ is multi- 
valued, one can only obtain to within a multiple of 2n, To 
obtain the actual frequency and motion corresponding to the 
root, p^, one sets Cj^=l and all other remaining C^=0 in Eqs. 

B-2) and (B-3) . Then, using the k eigenvector {B(0)}j^ from 
Eq. (B-7) as an initial condition, one would solve Eq. (B-1) with 

G=0 by numerical integration techniques for the resultant motion. 
Summarizing: To check for stability of a system of linear 

equations with periodic coefficients, obtain the eigenvalues 
of the "Transition Matrix" [Q] . If |Xj^|>l, one has instability. 
The traditional stability exponent Pj^ is related to through 
Eq?.(B-8) to (B-10). Two remarks on the above procedure should 
be noted. (1) The "Transition Matrix" [Q] can be formed by 
solving either the first order equations, Eqs. (B-1) with G=0, 
or the second order equations, Eqs. (A-1) with F=0 and periodic 
coefficients, whichever is more convenient for the integration 
scheme. (2) The above procedure will still apply even if the 
equations have constant coefficients. However, for such cases 
it is usually easier to form the matrix A given by Eq. (A-4) and 
obtain its eigenvalues Pj^ rather than to form the "Transition 
Matrix" [Q] and obtain its eigenvalues Aj^. 
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(b) Forced Response 

Solutions of Eq. (B-1 ) , or equivalently Eq. (A-1) with 
periodic coefficients, can be obtained by direct numerical 
integration using some convenient integration scheme. By 
proper choice of the initial conditions, one can eliminate 
all transients from the response and obtain the desired steady- 
state dynamic response by integrating through only one period T, 
instead of the very large number usually required to reach 
steady-state for lightly damped systems. A procedure for 
finding the proper initial conditions is given below. 

Solutions of Eq. (B-1) are of the general form, 

where homogeneous solution and i^p(t) is the 

particular solution. One can obtain a complete solution of 
Eq. (B-1) numerically for any given set of initial conditions. 

Call this solution y_g(t). One can add any number of additional 
homogeneous solutions Ay^j(t) having different initial conditions, 
to this solution. This would give a new solution to Eq. (B-1) , 

^(-t) = (B-12) 

which would have different initial conditions than those of 

Ze(^) . 

One can obtain all the homogeneous solutions of Eq. (B-1) 
by solving Eq. (B-1) with G=0 a total of 2N times, subject to 
the initial conditions y^”^ remaining y^=0, then Y 2 ~^ 

and all remaining y^=0, etc. In fact, this was done earlier 
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to investigate stability and resulted in the 2N homogeneous 
solutions (t) , , etc., respectively. Thus, one may 

write 



where [Q(t)] is the transition matrix at any instant of time, 
and C^, C 2 , ... are 2N arbitrary constants. The new solution 
Eq.(B-12) can be rewritten as 

-h LQC-t)]c (B-14) 

For a periodic solution over period T=2tt/^, one must have 
y^(T)=y(0). Placing Eq. (B-14) into this condition and solving 
for the arbitrary constants C gives, 

5ire(T) + [QWjc = y,E(o) + [Qto)]c 
c = [1 - MJ I 'ieir) - arsW] 

where it was noted that [Q(0)]=^, and [Q(T)]=[Q] is the 
"Transition Matrix" found earlier for the stability investigation. 
Placing these values of C back into Eq. (B-14) , the initial con- 
ditions for insuring a periodic solution become 

-I 

t [i - Q] { 


(B-16) 
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One can then solve the basic Eq. (B-1) numerically with these 
initial conditions to obtain a periodic solution over one period. 
It should be noted that if one had chosen the initial conditions 
for Xg(t) as would obtain simply, 

^Co) = [i " S] (B-17) 

This is particularly convenient form for finding the initial con- 
ditions for periodic solutions. 

An alternative form for determining the proper initial 
conditions for periodic solutions has been proposed by Friedmann 

and his coworkers^'^ in their work on wind turbines, namely, 

T 

-I f - -"I 

{^(Q> = [i - §] S. fW cit (B-18) 

This is similar to Eq. (B-17) , but does not use y . It seems 
easier to obtain ^.^(T) with initial conditions y^g(0)=0 and use 
Eq.(B-17), rather than obtaining [Q(t)] at every point and per- 
forming the indicated operations required by Eq. (B-18) .* 

The general procedure described by Eqs. (B-11) to (B-17) 
may be extended to deal also with nonlinear equations, 

i - Mt) ^ = r ("tj ) (B-19) 

where the right hand side now contains nonlinear functions of 
the coordinates. An iterative variation of the previous linear 
procedure to obtain the initial conditions for periodic solutions 

★ 

A comparison of these methods as well as a similar general de- 
rivation was given by Izadpanah (Ref. 13) . This was pointed out 
to the authors by Prof. D. A. Peters of Washington University, 
Saint Louis, Missouri. 
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of nonlinear equations is used by the MOSTAS Code^. The pro- 
cedure is as follows. First, a numerical solution is 

obtained to the nonlinear Eq. (B-19) for some estimate of the 
initial conditions Then each of the 2N elements of y.g(0) 

is perturbed a small amount and the resulting 2N solutions 
are obtained. This involves solving the nonlinear Eq. (B-19) 
subject to the initial conditions. 



which is in the same form as Eq. (B-13) . Then, again requiring 
the periodicity condition y^(T)=y^(0) and following through as 
before, will result in the same relation Eq. (B-16) found 
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previously. Because of the nonlinearities now present, the 
elements of [Q] as found from Eqs. (B-22) , (B-21) , (B-20) may 

vary with the amplitude of the initial condition used, Y.e ^ ^ i * 
This is in contrast to the linear case where [Q] remains always 
constant. Hence, an iterative application of Eq. (B-16) with a 
new corrected £g(0) should be done. If the nonlinearities are 
not too great, convergence to the required y.g(0) should be rapid. 

It should be remarked that the numerical procedure for 
forced response described in this section can also be used for 
the constant coefficient linear case, although it is probably 
easier there to obtain the solution by using Harmonic response 
methods given by Eqs. (A-6) to (A-10) . However, for cases where 
there is some nonlinearity, the present iterative approach 
becomes attractive. 

Summarizing, the Floquet methods described in this Appendix 
are based on a convenient numerical integration scheme and involve 
the computation of the "Transition Matrix," [Q] , from which both 
stability and the initial conditions for steady-state response 
solutions can be obtained. These methods seem attractive for 
large systems and can be modified to include nonlinearities in 
the equations. 
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Given a rotor with N blades rotating with rotation speed 

attached to a flexible tower. Because the tower motions x^ 

are described in a fixed reference frame while the blade motions 

are described relative to a rotating frame, the resulting 

equations may have mass, damping, or stiffness coefficients 

th 

which are functions of the azimuthal position of the k blade 
A typical such set of equations is given, for example, in 
Refs. 4 and 14 as, 

2 . 

Mx ^ ^ 7. I == F^(^) 

d't 

(C-1) 

( ^ = 1 > "2.^ N ) 

where the azimuthal position is, 

^ = 5lt -t (.^-|)27 t/n 

The first equation above represents force equilibrium for the 

tower motion x, while the remaining N equations represent force 

(k) 

equilibrium for the motion of each of the N blades B . The 

above equations are readily generalized to more tower motions 

( k ) 

x^, and more blade coordinates for each blade B^ 
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(a) Stability 

To examine Eqs.(C-l) for stability, one sets F =0 and 
( k ) 

Fg =0 to obtain homogeneous equations. 

For rotors with 3 or more blades N_>3, one may eliminate 
the periodic coefficients in these equations by introducing new 


multiblade coordinates b^( t) , . b^^ ( t) , b^^ ( t) , b 2 g ( t) , ... 
that 


such 





where the total number of coordinates b^ taken is equal to the 

number of blades, N. Note, if N = even number, then the 

k ““ X 

additional mode (-1) b^ is needed to complete the set. The 

above form is an equivalent modal representation of the N blades. 

Substituting these into Eqs. (C-1) , then multiplying the last N 

k-1 

equations by sin cos 1, sin cos ••• (“D 

respectively, then slimming these last N equations and noting the 
trigonometric identities. 


=. Gya- ^ ~ 'z. + ^ 

(Gr<2. 4;^ + 

C^r8- ^ C£><z.(iM + h) 
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Nl 


^ 1 ^ V|^ 

A - 1 


k; 

Jk.= 

4-1 




L (-0 - 


^=1 

N 


Z (' iZZii-wv = 





+ 0 ^ Wl 

= ZNj .. - 

0 

■f-or >T^ 

Kl, 2 .N,--- 

N 0 ?Q. 

-f-o r 

= N ^ XN^ ... 

0 

4 -or no 

N j ZN^--' 

r i\i v*\ 4 ^ , 


101 =. N/i. ^ 3 h ^2 ^ ' 

1 0 

•for 

m .if lV/ 2 . ^ bN/z y . 

r N ynti 

for 

iv\ = N/Z j 3 S/l^ 

1 0 

for 

wo N/z j 3 M/ 2 ^ 


(C- 4 ) 


results in a new set of differential equations in the variables 

X , ^Ic' ^ 2 s ' ^ 2 c *** ^A' ^ blades, now all 

have constant coefficients namely, 


Mr + C^r ' b ,^ = o 

1 b„ Cjb„ + 25211;,^- 51C^b„ 

+ 2su(;„ + sic,b,, + 1 b,,-^ 


N 

2 - 


= 0 


ii 

2 - 

N 


= 0 


[t bV-t Cp tv 


= 0 


il 

2 - 


J b^5 -h Cg + i'^'~ °X5 ^*2c"~ ^2C 


= 0 


N 




(C-5) 


0 
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These equations may then be investigated for stability using 

the standard constant coefficient techniques described earlier. 

It should be noted that because of the form of Eqs.(C-l), the 

equations above uncouple into several smaller groups, and in 

fact only the first three equations above are coupled together 

and involve the tower motion variable x. For additional details 

and applications of multiblade coordinates, see Hohenemser and 
3 15 

Yin and Johnson . Multiblade coordinates were originally 

1 6 

introduced by Coleman and Feingold in their studies of heli- 
copter ground resonance. 

For rotors with 2 blades, N=2, the analysis is more 
difficult because the rotor disk no longer has polar symmetry. 
For this case, one may use multiblade coordinates together with 
harmonic balance methods to arrive at approximate solutions. 

The multiblade transformation Eq. (C-3) for 2 blades can be ex- 
pressed in terms of coordinates b^(t) and b^(t) as, 

(C-6) 

These coordinates are introduced into Eqs.(C-l), then summing 
and subtracting the last two equations of Eqs.(C-l) respectively 
while noting that sin = -sin >1;^^ and cos 1^2 “ -cos results 
in a new set of differential equations in the variables x, b^, 
b^ v/hich still have periodic coefficients. Then, expanding each 
of the coordinates in a harmonic series. 



-53 


^ = "X*, + -I- 00 ^ Sit -h -X^5 XUvu 2 Sl-t + " • 

(C-7) 

t? = b b Ax^^l't- + b_ C^52t + b_ 2 SZt + - - • 

T ‘^TO TIS TIC I 2S 

^AIS ^ ^-is 

where x^, ^ 2 . 3 ' ^10' ^TlS' *'* functions of time, then 

placing these into the previous equations and balancing out 
each haimonic term in each equation, will yield an infinite set 
of constant coefficient differential equations which can be 
truncated at some point for approximate solutions. These 
truncated equations may again be examined for stability using 
the standard constant coefficient techniques described earlier. 

Often, depending on the form of Eqs.(C-l), the resulting 
constant coefficient differential equations will uncouple into 
several smaller coupled systems of equations which may be examined 
independently of one another. For example, for the case of 
Eqs. (C-1) , one smaller coupled system would involve the variables 

^0' ^2S' ^2C' ^AlS' ^AlC' •*' another would involve x^^, 

^IS' ^AO' ^A2S' '^A2C' *'* such systems, one could use an 

alternate extended form of the multiblade coordinate transformation 
Eq. (C-3) namely, 

'X = -«UAA-Z51t f 'X^^cer>~2S± + 

(C-8) 
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together with the harmonic balance method to solve the problem. 

This works here, since the form given by Eq. (C-8) exactly 
duplicates the motion of the two blades given by the general 
case Eqs. (C-6) and (C-7) , since sin ~ “Sin cos 4^2 “ 

-cos 4^^, and only Xq, ^ ^AlC' **’ present. 

4 

See Sheu for an application of this alternate extended form of 
the multiblade transformation Eq. (C-8) to a simple two-bladed 
rotor in ground resonance, Eqs. (C-1) . Solutions involving as 
little as three terms, x^, gave reasonable approximations 

to the primary instability regions. However, in more general 
cases (for example, if the first equation of Eqs. (C-1) had an 
additional terms M^x cos 4'^_ or k^x cos 4^^ present) , the resulting 
equations would not split into two smaller groups, and the general 
harmonic balance method Eqs. (C-6) and (C-7) would have to be used. 

Indeed, for the more general case mentioned above, one 
would also investigate the system for direct Mathieu equation 
type instabilities of half integer order f2/2, 3f2/2, ... by in- 
troducing additional harmonic terms sin mf?t and cos mPt where 
m=l/2, 3/2, 5/2, ... into Eqs. (C-7), and harmonically balancing 
as before. These terms would not couple in with the previous 
equations and can be solved independently of them. The primary 

17 

instability region would result from the n/2 terms. See Bolotin 
for further details of the general harmonic balance method. 
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"K — ocro^ (O^t - 

* ( SsR ~ ^SI'^M,'^'} (C-11) 

■" ( •=’icR - t|cr-^‘<^>"'<-) 

^ ( ^ZSR ~ ^2SI 4- . . , 


The tower thus oscillates at frequency u in the fixed frame 

m 

whereas the blades may oscillate at frequencies w , w + n, 

' m Tti 

u - n, u + 2n, u - 2fi, ... relative to the rotating frame, 

depending on which coordinates are excited. For example, 

in the case of a 3 bladed rotor N=3, the term of Eq. (C-10) 

leads to a constant term on the right-hand-side which excites 

the X, b^^ coordinates and results in a constant x and a 

(k) 

3 frequency of while the leads to a cos 3i(; term 

on the right-hand-side which again excites x, b^^ and 

(k) 

results in an x frequency of 3Q and 6 frequencies of 29 ,, 


For rotors with 2 blades N=2, one can use the harmonic 
balance methods of the previous section. The steady-state 
periodic tower and blade forces given by Eqs. (C-10) are sub- 
stituted into the basic equations of motions [Eqs.(C-l)]. 


One then introduces the new coordinates given by Eqs. (C-6 ) , 
then sums and subtracts the last two blade equations, then 
expands the tower and blade motions as given by Eq. (C-7) , only 
now the coordinates Xq, ^IS^ ^TO^ ^TlS^ • • • 


etc. are taken 
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to be constants rather than functions of time. Harmonically 
balancing the various terms in each equation results in a 
truncated set of algebraic equations which can be solved to 
obtain the coordinates ^Is' ^To' **' corresponding 

to the given forcing excitations ^^Is' ^Bo' ^Bo' ^Bls' 

etc. The resulting tower and blade motions are then given 
directly by Eqs. (C-7) and (C-6) . The resulting set of algebraic 
equations will often uncouple into smaller coupled sets of 
equations which can be examined independently of one another. 

This procedure is similar to that for the constant coefficient 
forced response case Eq. (A-8) , except now, the periodic co- 
efficients couple the different harmonics together. Thus, the 
solution will consist of many harmonics mf^ even if only one 
forcing harmonic were present alone. 

Summarizing, the multiblade coordinate and harmonic balance 
methods described in this Appendix involve first the introduction 
of multiblade coordinates in order to take out the periodic co- 
efficients from the blades [Eqs.(C-3) for N > 3 ] , or to obtain 
a better ordered system of equations [Eqs. (C-6) for N=2]. Then 
harmonic balance methods Eqs. (C-7) are used to deal with any 
remaining periodic coefficients. These methods seem attractive 
for smaller systems and can give considerable insight into the 
origin and nature of instabilities and the various harmonics 
present in the forced response. 
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Appendix D 
Rotating Coordinates 


As an addendum to the previous multiblade coordinates and 
harmonic balance methods, it should be mentioned that for some 
problems, the use of rotating coordinates is also convenient. 
For example, in the case of a 2-bladed rotor on isotropic tower 
supports (same tower mass, damping, and stiffness in two direc- 
tions, x^ and ^ 2 ) , Eqs.(C-l) would read. 


Mi', t qa, 




,l -fe- 


(D-1) 




JV 




N) 


One can then express the tower motions x^ and X 2 in terms of 
rotating coordinates ^2 rotate with the rotor, as 

•X| — c^rcL. xwu- 52r 

“Xj, = “ Pi rifrd. 52 . 't (D-2) 

where the rotation = f^t is taken from the X 2 axis towards 
the axis. Placing these equations into Eqs.(D-l) , then 
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multiplying the first two equations by cos and sin 
respectively and subtracting, then multiplying the first two 
equations by sin and cos and adding, then subtracting 
the third and fourth equations, then adding the third and 
fourth equations will result in a new set of differential 
equations in the variables ^ 2 ' which now all have 

constant coefficients, namely, 

+ I 52t - S2t 

-4-^siL'^ = F^I Ccr^ Sl-t 


(i> 2It> 2C^ F^ 


(') 

3 


ZtL], + = f /'^+ 


In the above, b^ = and b^ = [3^^^ - 3^^h/2 

are the same coordinates introduced earlier in Eqs.(C-6). These 
differential equations may then be investigated for stability 
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and forced response using the standard constant coefficient 

techniques described earlier. Such analyses of a 2-bladed 

rotor on isotropic tower supports were also performed by Coleman 
16 

and Feingold in their studies of helicopter ground resonance. 

Rotating coordinates are often used in rotating machinery 

shaft critical speed problems, and are useful for dealing with 

problems of rotors with unsymmetrical mass, unsymmetrical damping, 

or unsymmetrical shaft stiffness supported on isotropic bearings. 

1 8 

See for example, Bolotin . For such problems, one can readily 
set up the equations of motion in the rotating frame directions, 
and the fixed supports will introduce no periodic terms because 
of their isotropic nature. For vertical axis wind turbines, 
such rotating coordinates for the blades are useful since the 
tower supports are generally isotropic due to the symmetrically 
arranged guy wires. For horizontal axis wind turbines, the 
tower supports are generally not isotropic, hence periodic co- 
efficients will remain in the equations when using rotating 
coordinates. If the support anisotropy is not too large, one 
can again additionally introduce harmonic balance methods to 
eliminate the periodic coefficients, as was done in the previous 
section. 

Summarizing, the rotating coordinates described in this 
Appendix can also be used to effectively eliminate the periodic 
coefficients in problems involving unsymmetrical rotors on 
isotropic tower supports. These can often be used in rotating 
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shaft critical speed problems and for vertical axis wind 
turbines. If the support anisotropy is not too large, harmonic 
balance methods may additionally be used to deal with any 
remaining periodic coefficients. 
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INPUTS 

Steady Wind 
Wind Shear 
Tov;er Shadow 
Cl-ravity 
Yav; 

Gust 

Stat. Unbal. 




f — Total rotor load at fixed shaft 
o 

f — Perturbation pod load applied to rotor 

f, — Fixed shaft blade load R, — Fixed shaft blade coord, 

b 

f -- Perturbation blade load x — Rotor shaft displace. 

f„^„ — f, + f — Rotation speed 

TOT b ^ 


FIG. 1 


Layout of MOSTAS Computer Code 
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FIG. 2 


WINDLASS System Block Diagram 
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